/************************************************************************
 * MechSys - Open Library for Mechanical Systems                        *
 * Copyright (C) 2013 Sergio Torres                                     *
 *                                                                      *
 * This program is free software: you can redistribute it and/or modify *
 * it under the terms of the GNU General Public License as published by *
 * the Free Software Foundation, either version 3 of the License, or    *
 * any later version.                                                   *
 *                                                                      *
 * This program is distributed in the hope that it will be useful,      *
 * but WITHOUT ANY WARRANTY; without even the implied warranty of       *
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the         *
 * GNU General Public License for more details.                         *
 *                                                                      *
 * You should have received a copy of the GNU General Public License    *
 * along with this program. If not, see <http://www.gnu.org/licenses/>  *
 ************************************************************************/

// STL
#include <iostream>
#include <cmath>

// MechSys
#include <mechsys/dem/domain.h>
#include <mechsys/util/fatal.h>
#include <mechsys/util/array.h>
#include <mechsys/util/util.h>
#include <mechsys/linalg/matvec.h>
#include <mechsys/mesh/structured.h>
#include <mechsys/mesh/unstructured.h>

using std::cout;
using std::endl;
using std::ifstream;
using std::min;
using std::max;


struct UserData
{
    double sx;
    double sy;
    double p0;
    double  p;
    double  Q;         
    double Bn;
    double Ka;         
    double Kb;         
    double V0;
    double V0f;
    double Vf;
    double ev;
    double To;
    double ev0;
    double Vref;
    DEM::Particle       * px0;
    DEM::Particle       * py0;
    DEM::Particle       * pz0;
    DEM::Particle       * px1;
    DEM::Particle       * py1;
    DEM::Particle       * pz1;
    std::ofstream      oss_ss;       ///< file for stress strain data
};

void Setup1 (DEM::Domain & dom, void * UD)
{
    UserData & dat = (*static_cast<UserData *>(UD));
    double areax = (dat.pz1->x(2) - dat.pz0->x(2))*(dat.py1->x(1) - dat.py0->x(1));
    double areay = (dat.pz1->x(2) - dat.pz0->x(2))*(dat.px1->x(0) - dat.px0->x(0));
    dat.px0->Ff  =  areax*dat.sx*OrthoSys::e0;
    dat.px1->Ff  = -areax*dat.sx*OrthoSys::e0;
    dat.py0->Ff  =  areay*dat.sy*OrthoSys::e1;
    dat.py1->Ff  = -areay*dat.sy*OrthoSys::e1;
    
    dat.px0->Ff  =  dat.px0->Ff*min(1.0,dom.Time/50.0); 
    dat.px1->Ff  =  dat.px1->Ff*min(1.0,dom.Time/50.0); 
    dat.py0->Ff  =  dat.py0->Ff*min(1.0,dom.Time/50.0);
    dat.py1->Ff  =  dat.py1->Ff*min(1.0,dom.Time/50.0);
    
     double volumecontainer = (dat.px1->x(0)-dat.px0->x(0)-dat.px1->Props.R+dat.px0->Props.R)*
                             (dat.py1->x(1)-dat.py0->x(1)-dat.py1->Props.R+dat.py0->Props.R)*
                             (dat.pz1->x(2)-dat.pz0->x(2)-dat.pz1->Props.R+dat.pz0->Props.R);

    // Here is the formula for the pressure as a function of ev
    dat.Vf  = volumecontainer - dat.V0 + dat.V0f;
    dat.ev  = (dat.Q*dom.Time + dat.V0f - dat.Vf)/(dat.V0f);
    dat.p   = dat.p0*min(1.0,dom.Time/50.0);
    
     for (size_t i=0;i<dom.BInteractons.Size();i++)
    {
        DEM::BInteracton * BI = dom.BInteractons[i];
        if (!BI->valid) 
        {
            //std::cout << "0" << " " << BI->P1->Tag << " " << BI->P2->Tag << std::endl;

            BI->P1->Tag = 1;
            BI->P2->Tag = 1;

            //std::cout << "1" << " " << i << std::endl;

            Vec3_t n1,n2;
            double a1,a2;
            BI->P1->Faces[BI->IF1]->Normal(n1);
            BI->P2->Faces[BI->IF2]->Normal(n2);
            a1 = BI->P1->Faces[BI->IF1]->Area();
            a2 = BI->P2->Faces[BI->IF2]->Area();

            //std::cout << "2" << " " << BI->P1->Tag << " " << BI->P2->Tag << std::endl;

            Vec3_t F1 = -dat.p*a1*n1;
            Vec3_t F2 = -dat.p*a2*n2;
            //std::cout << "2" << " " << F1 << " " << F2 << std::endl;

            BI->P1->F += F1;
            BI->P2->F += F2;

            Vec3_t c1,c2;
            BI->P1->Faces[BI->IF1]->Centroid(c1);
            BI->P2->Faces[BI->IF2]->Centroid(c2);

            Vec3_t x1,x2;
            x1 = c1 - BI->P1->x;
            x2 = c2 - BI->P2->x;

            Vec3_t T, Tt;
            Tt = cross (x1,F1);
            Quaternion_t q1;
            Conjugate (BI->P1->Q,q1);
            Rotation  (Tt,q1,T);
            BI->P1->T -= T;
            Tt = cross (x2,F2);
            Quaternion_t q2;
            Conjugate (BI->P2->Q,q2);
            Rotation  (Tt,q2,T);
            BI->P2->T += T;
         }
    }
}

void Setup2 (DEM::Domain & dom, void * UD)
{
    UserData & dat = (*static_cast<UserData *>(UD));
    double areax = (dat.pz1->x(2) - dat.pz0->x(2))*(dat.py1->x(1) - dat.py0->x(1));
    double areay = (dat.pz1->x(2) - dat.pz0->x(2))*(dat.px1->x(0) - dat.px0->x(0));
    dat.px0->Ff  =  areax*dat.sx*OrthoSys::e0;
    dat.px1->Ff  = -areax*dat.sx*OrthoSys::e0;
    dat.py0->Ff  =  areay*dat.sy*OrthoSys::e1;
    dat.py1->Ff  = -areay*dat.sy*OrthoSys::e1;

    //double volumecontainer = (dat.px1->x(0)-dat.px0->x(0)-dat.px1->Props.R+dat.px0->Props.R)*
                             //(dat.py1->x(1)-dat.py0->x(1)-dat.py1->Props.R+dat.py0->Props.R)*
                             //(dat.pz1->x(2)-dat.pz0->x(2)-dat.pz1->Props.R+dat.pz0->Props.R);
    double volumefracture = 0.0;
    for (size_t i=0;i<dom.BInteractons.Size();i++) 
    {
        if (!dom.BInteractons[i]->valid)
        {
            //volumefracture += norm(dom.BInteractons[i]->Fnet)*dom.BInteractons[i]->L0/dat.Bn;
            DEM::Particle * P1 = dom.BInteractons[i]->P1;
            DEM::Particle * P2 = dom.BInteractons[i]->P2;
            Vec3_t c1,c2;
            P1->Faces[dom.BInteractons[i]->IF1]->Centroid(c1);
            P2->Faces[dom.BInteractons[i]->IF2]->Centroid(c2);
            double delta = norm(c2-c1)-P1->Props.R-P2->Props.R;
            if (delta<0.0) delta = 0.0;
            volumefracture += delta*dom.BInteractons[i]->Area*(1.0+4.0*P1->Props.R/sqrt(dom.BInteractons[i]->Area));
        }
    }

    // Here is the formula for the pressure as a function of ev
    //dat.Vf  = volumecontainer - dat.V0 + dat.V0f;
    dat.Vf  = volumefracture + dat.V0f;
    dat.ev  = (dat.Q*(dom.Time-dat.To) + dat.V0f - dat.Vf)/(dat.Q*(dom.Time-dat.To) + dat.V0f);
    dat.p   = (dat.p0 + dat.Kb*dat.ev)/(1.0-dat.Ka*dat.ev);
    dat.p   = max(dat.p,0.0);

    //In this point you change the formula for the pressure

#ifdef USE_OMP
    #pragma omp parallel for schedule (static) num_threads(dom.Nproc)
#endif
    for (size_t i=0;i<dom.BInteractons.Size();i++)
    {
        DEM::BInteracton * BI = dom.BInteractons[i];
        if (!BI->valid) 
        {
            //std::cout << "0" << " " << BI->P1->Tag << " " << BI->P2->Tag << std::endl;

            BI->P1->Tag = 1;
            BI->P2->Tag = 1;

            //std::cout << "1" << " " << i << std::endl;

            Vec3_t n1,n2;
            double a1,a2;
            BI->P1->Faces[BI->IF1]->Normal(n1);
            BI->P2->Faces[BI->IF2]->Normal(n2);
            a1 = BI->P1->Faces[BI->IF1]->Area();
            a2 = BI->P2->Faces[BI->IF2]->Area();

            //std::cout << "2" << " " << BI->P1->Tag << " " << BI->P2->Tag << std::endl;

            Vec3_t F1 = -dat.p*a1*n1;
            Vec3_t F2 = -dat.p*a2*n2;
            //std::cout << "2" << " " << F1 << " " << F2 << std::endl;


            Vec3_t c1,c2;
            BI->P1->Faces[BI->IF1]->Centroid(c1);
            BI->P2->Faces[BI->IF2]->Centroid(c2);

            Vec3_t x1,x2;
            x1 = c1 - BI->P1->x;
            x2 = c2 - BI->P2->x;

            Vec3_t T, Tt;
            Tt = cross (x1,F1);
            Quaternion_t q1;
            Conjugate (BI->P1->Q,q1);
            Rotation  (Tt,q1,T);
            omp_set_lock  (&dom.BInteractons[i]->P1->lck);
            BI->P1->F += F1;
            BI->P1->T -= T;
            omp_unset_lock(&dom.BInteractons[i]->P1->lck);
            Tt = cross (x2,F2);
            Quaternion_t q2;
            Conjugate (BI->P2->Q,q2);
            Rotation  (Tt,q2,T);
            omp_set_lock  (&dom.BInteractons[i]->P2->lck);
            BI->P2->F += F2;
            BI->P2->T += T;
            omp_unset_lock(&dom.BInteractons[i]->P2->lck);

            //std::cout << "3" << std::endl;
            //Vec3_t Xmin,Xmax;
            //dom.BoundingBox(Xmin,Xmax);
            //Vec3_t DomCtr = 0.5*(Xmin+Xmax);
            //std::cout << "1 " << BI->P1->x - DomCtr << " " << F1 << std::endl;
            //std::cout << "2 " << BI->P2->x - DomCtr << " " << F2 << std::endl;
        }
    }
}

void Report (DEM::Domain & dom, void *UD)
{
    UserData & dat = (*static_cast<UserData *>(UD));
    if (dom.idx_out==0)
    {
        String fs;
        fs.Printf("%s_walls.res",dom.FileKey.CStr());
        dat.oss_ss.open(fs.CStr());
        dat.oss_ss << Util::_10_6 << "Time" << Util::_8s << "Volfrac" << Util::_8s << "ev" << Util::_8s << "Pressure" << Util::_8s << "Displc_x" << Util::_8s << "Displc_y"  << Util::_8s << "p_x" << Util::_8s << "p_y" << std::endl;
    }
    if (!dom.Finished) 
    {
        String ff;
        ff.Printf     ("%s_frac_%04d",dom.FileKey.CStr(), dom.idx_out);
        dom.WriteFrac (ff.CStr());
        ff.Printf     ("%s_bf_%04d",dom.FileKey.CStr(), dom.idx_out);
        dom.WriteBF   (ff.CStr());
		double areax = (dat.pz1->x(2) - dat.pz0->x(2))*(dat.py1->x(1) - dat.py0->x(1));
        double areay = (dat.pz1->x(2) - dat.pz0->x(2))*(dat.px1->x(0) - dat.px0->x(0));
        dat.oss_ss << Util::_10_6 << dom.Time << Util::_8s << dat.Vf << Util::_8s << dat.ev << Util::_8s << dat.p << Util::_8s << dat.px1->x(0)-dat.px0->x(0) << Util::_8s << dat.py1->x(1)-dat.py0->x(1) << Util::_8s << dat.px0->Ff(0)/areax << Util::_8s << dat.py0->Ff(1)/areay << std::endl;
    }
    else
    {
        dat.oss_ss.close();
    }
}

void GenFromMesh(DEM::Domain & dom, char const * VertsFile, char const * CellsFile, double thickness, double scale, double R, double rho, bool Cohesion = false)
{
    String fnv(VertsFile); fnv.append(".mesh");
    String fnf(CellsFile); fnf.append(".mesh");
    ifstream fnvf(fnv.CStr());
    ifstream fnff(fnf.CStr());
    if (!Util::FileExists(fnv)) throw new Fatal("File <%s> not found",fnv.CStr());
    if (!Util::FileExists(fnf)) throw new Fatal("File <%s> not found",fnf.CStr());
    size_t ncol = 0;
    Vec3_t Vtemp = OrthoSys::O;
    size_t number;
    Array<Vec3_t> Vtotal;
    while (!fnvf.eof())
    {
        if (ncol==0) fnvf >> number;
        if (ncol==1) fnvf >> Vtemp(ncol-1);
        if (ncol==2) fnvf >> Vtemp(ncol-1);
        ncol++;
        ncol = ncol%3;
        if (ncol==0)
        {
            Vtotal.Push(scale*Vtemp);
            //std::cout << Vtotal[Vtotal.Size()-1] << std::endl;
        }
    }
    Array<Array <size_t> > Ctotal;
    Array<size_t>        Cpartial;
    size_t ncoltot = 0;
    ncol = 0;
    while (!fnff.eof())
    {
        if (ncol==0)
        {
            fnff >> ncoltot;
            ncoltot++;
            Cpartial.Resize(0);
            Cpartial.Resize(ncoltot-1);
        }
        else
        {
            fnff >> number;
            Cpartial[ncol-1] = number;
        }
        ncol++;
        ncol = ncol%ncoltot;
        if (ncol==0)
        {
            Ctotal.Push(Cpartial);
            //std::cout << Ctotal[Ctotal.Size()-1][0] << " " << Ctotal[Ctotal.Size()-1][1] << " " << Ctotal[Ctotal.Size()-1][2] << " " << Ctotal[Ctotal.Size()-1][3] << " " << std::endl;
        }
    }

    //for (size_t i=0;i<Ctotal.Size();i++)
    //{
        //for (size_t j=0;j<Ctotal[i].Size();j++)
        //{
            //std::cout << Ctotal[i][j] << " ";
        //}
        //std::cout << std::endl;
    //}

    //return;


    size_t IIndex = dom.Particles.Size();
    Array<Vec3_t> V;
    Array<Array <int> > E;
    Array<Array <int> > F;
    for (size_t i=0;i<Ctotal.Size();i++)
    {
        size_t nverts = Ctotal[i].Size();
        V.Resize(2*nverts);
        E.Resize(3*nverts);
        F.Resize(nverts+2);
        //std::cout << "1" << std::endl;
        for (size_t j=0;j<Ctotal[i].Size();j++)
        {
            V[j+nverts] = Vtotal[Ctotal[i][j]] + Vec3_t(0.0,0.0,-thickness);
            V[j]        = Vtotal[Ctotal[i][j]];
        }
        //std::cout << "2" << std::endl;
        for (size_t j=0;j<nverts;j++)
        {
            E[j].Push(j);
            E[j].Push((j+1)%nverts);
            E[j+nverts].Push(j);
            E[j+nverts].Push(j+nverts);
            E[j+2*nverts].Push(j+nverts);
            E[j+2*nverts].Push((j+1)%nverts+nverts);
        }
        //std::cout << "3" << std::endl;
        for (size_t j=0;j<nverts;j++)
        {
            F[j].Push(j);
            F[j].Push(j+nverts);
            F[j].Push((j+1)%nverts+nverts);
            F[j].Push((j+1)%nverts);
            F[nverts].Push(j);
            F[nverts+1].Push(2*nverts-j-1);
        }
        double vol; // volume of the polyhedron
        Vec3_t CM;  // Center of mass of the polyhedron
        Mat3_t It;  // Inertia tensor of the polyhedron
        DEM::PolyhedraMP(V,F,vol,CM,It);
        DEM::Erosion(V,E,F,R);
        dom.Particles.Push (new DEM::Particle(-1, V,E,F,OrthoSys::O,OrthoSys::O,R,rho));
        dom.Particles[dom.Particles.Size()-1]->Eroded = true;
        dom.Particles[dom.Particles.Size()-1]->x       = CM;
        dom.Particles[dom.Particles.Size()-1]->Props.V = vol;
        dom.Particles[dom.Particles.Size()-1]->Props.m = vol*rho;
        Vec3_t I;
        Quaternion_t Q;
        Vec3_t xp,yp,zp;
        Eig(It,I,xp,yp,zp);
        CheckDestroGiro(xp,yp,zp);
        I *= rho;
        Q(0) = 0.5*sqrt(1+xp(0)+yp(1)+zp(2));
        Q(1) = (yp(2)-zp(1))/(4*Q(0));
        Q(2) = (zp(0)-xp(2))/(4*Q(0));
        Q(3) = (xp(1)-yp(0))/(4*Q(0));
        Q = Q/norm(Q);
        dom.Particles[dom.Particles.Size()-1]->I     = I;
        dom.Particles[dom.Particles.Size()-1]->Q     = Q;
        double Dmax = DEM::Distance(CM,V[0])+R;
        for (size_t i=1; i<V.Size(); ++i)
        {
            if (DEM::Distance(CM,V[i])+R > Dmax) Dmax = DEM::Distance(CM,V[i])+R;
        }
        dom.Particles[dom.Particles.Size()-1]->Ekin = 0.0;
        dom.Particles[dom.Particles.Size()-1]->Erot = 0.0;
        dom.Particles[dom.Particles.Size()-1]->Dmax  = Dmax;
        dom.Particles[dom.Particles.Size()-1]->PropsReady = true;
        dom.Particles[dom.Particles.Size()-1]->Index = dom.Particles.Size()-1;
    }

    if(Cohesion)
    {
        double tol1 = 1.0e-8;
        double tol2 = 1.0e-3;
        for (size_t i=IIndex;i<dom.Particles.Size()-1;i++)
        {
            DEM::Particle * P1 = dom.Particles[i];
            for (size_t j=i+1;j<dom.Particles.Size();j++)
            {
                DEM::Particle * P2 = dom.Particles[j];
                if (DEM::Distance(P1->x,P2->x)<P1->Dmax+P2->Dmax)
                {
                    for (size_t k=0;k<P1->Faces.Size();k++)
                    {
                        DEM::Face * F1 = P1->Faces[k];
                        Vec3_t n1,c1;
                        F1->Normal  (n1);
                        F1->Centroid(c1);
                        bool found = false;
                        for (size_t l=0;l<P2->Faces.Size();l++)
                        {
                            DEM::Face * F2 = P2->Faces[l];
                            Vec3_t n2,c2;
                            F2->Normal  (n2);
                            F2->Centroid(c2);
                            Vec3_t n = 0.5*(n1-n2);
                            n/=norm(n);
                            if ((fabs(dot(n1,n2)+1.0)<tol1)
                               &&(fabs(Distance(c1,*F2)-2*R)<tol2)
                               &&(fabs(Distance(c2,*F1)-2*R)<tol2))
                            {
                                dom.BInteractons.Push(new DEM::BInteracton(P1,P2,k,l));
                                found = true;
                                break;
                            }
                        }
                        if (found) break;
                    }
                }
            }
        }
        
    }
}

int main(int argc, char **argv) try
{
    
    if (argc<2) throw new Fatal("This program must be called with one argument: the name of the data input file without the '.inp' suffix.\nExample:\t %s filekey\n",argv[0]);

    UserData dat;
    DEM::Domain dom(&dat);
    dom.Dilate = true;
    size_t Nproc = 1; 
    if (argc==3) Nproc=atoi(argv[2]);
    String filekey  (argv[1]);
    String filename (filekey+".inp");
    if (!Util::FileExists(filename)) throw new Fatal("File <%s> not found",filename.CStr());
    ifstream infile(filename.CStr());


    String cellsfile;   // File containing the cell information
    String vertsfile;   // File containing the vertices information
    double scale;       // Scale the mesh particles    
    double Kn;          // Normal stiffness
    double Kt;          // Tangential stiffness
    double Gn;          // Normal dissipative coefficient
    double Gt;          // Tangential dissipative coefficient
    double Mu;          // Microscopic friction coefficient
    double Bn;          // Cohesion normal stiffness
    double Bt;          // Cohesion tangential stiffness
    double Bm;          // Cohesion torque stiffness
    double Eps;         // Cohesion torque stiffness
    double R;           // Spheroradius
    double sx;          // Stress in the x direction
    double sy;          // Stress in the y direction
    double p0;          // Fluid initial pressure
    double Q;           // Flow rate in Vol/time
    double Ka;          // Slope of the pressure coeffcient
    double Kb;          // Free constant of the pressure coefficient.
    double dt;          // Time step
    double dtOut;       // Time step for output
    double To;          // Time for the hydrostatic test
    double Tf;          // Time for the test

    infile >> cellsfile;  infile.ignore(200,'\n');
    infile >> vertsfile;  infile.ignore(200,'\n');
    infile >> scale;      infile.ignore(200,'\n');
    infile >> Kn;         infile.ignore(200,'\n');
    infile >> Kt;         infile.ignore(200,'\n');
    infile >> Gn;         infile.ignore(200,'\n');
    infile >> Gt;         infile.ignore(200,'\n');
    infile >> Mu;         infile.ignore(200,'\n');
    infile >> Bn;         infile.ignore(200,'\n');
    infile >> Bt;         infile.ignore(200,'\n');
    infile >> Bm;         infile.ignore(200,'\n');
    infile >> Eps;        infile.ignore(200,'\n');
    infile >> R;          infile.ignore(200,'\n');
    infile >> sx;         infile.ignore(200,'\n');
    infile >> sy;         infile.ignore(200,'\n');
    infile >> p0;         infile.ignore(200,'\n');
    infile >> Q;          infile.ignore(200,'\n');
    infile >> Ka;         infile.ignore(200,'\n');
    infile >> Kb;         infile.ignore(200,'\n');
    infile >> dt;         infile.ignore(200,'\n');
    infile >> dtOut;      infile.ignore(200,'\n');
    infile >> To;         infile.ignore(200,'\n');
    infile >> Tf;         infile.ignore(200,'\n');

    dom.Alpha = 0.5*R;
    if (cellsfile=="voronoi")
    {
        dom.AddVoroPack(-1, R, scale, scale, 1.0, size_t(scale), size_t(scale), 1, 3000.0, true, false, 1000, 1.0, Vec3_t(0.0,0.0,1.0));
    }
    else
    {
        GenFromMesh(dom,vertsfile.CStr(),cellsfile.CStr(),1.0,scale,R,3000.0,true);
    }


    dom.GenBoundingBox (/*InitialTag*/-2, 0.2*R, /*Cf*/1.3,false);
    
    dat.sx  = sx;
    dat.sy  = sy;
    dat.p0  = p0;
    dat.p   = p0;
    dat.Bn  = Bn;
    dat.Q   = Q; 
    dat.Ka  = Ka;
    dat.Kb  = Kb;
    dat.To  = To;

    dat.px0 = dom.GetParticle(-3);
    dat.px1 = dom.GetParticle(-2);
    dat.py0 = dom.GetParticle(-5);
    dat.py1 = dom.GetParticle(-4);
    dat.pz0 = dom.GetParticle(-7);
    dat.pz1 = dom.GetParticle(-6);
    dat.px0->FixVeloc();
    dat.px1->FixVeloc();
    dat.py0->FixVeloc();
    dat.py1->FixVeloc();
    dat.pz0->FixVeloc();
    dat.pz1->FixVeloc();
    dat.px0->vxf = false;
    dat.px1->vxf = false;
    dat.py0->vyf = false;
    dat.py1->vyf = false;

    Vec3_t Xmin,Xmax;
    dom.BoundingBox(Xmin,Xmax);
    Vec3_t DomCtr = 0.5*(Xmin+Xmax);
    
    double dist = norm(Xmin-Xmax);
    DEM::Particle * ParCtr = dom.Particles[0];
    for(size_t i=0;i<dom.Particles.Size();i++)
    {
        if ((DEM::Distance(DomCtr,dom.Particles[i]->x) < dist)&&(dom.Particles[i]->IsFree()))
        {
            dist   = DEM::Distance(DomCtr,dom.Particles[i]->x);
            ParCtr = dom.Particles[i];
        }
    }
    ParCtr->Tag = 1000;
    Array<int> DeleteTags(1);
    DeleteTags[0]  = 1000;
    dom.DelParticles(DeleteTags);

    for(size_t i=0;i<dom.Particles.Size();i++) dom.Particles[i]->Index = i;
    
    double Bno = Bn;          // Cohesion normal stiffness
    double Bto = Bt;          // Cohesion tangential stiffness
    
    Dict B;
    B.Set(-1,"Kn Kt Gn Gt Mu Bn Bt Bm Eps",Kn,Kt,Gn,Gt,Mu ,Bno,Bto ,Bm ,    Eps);
    B.Set(-2,"Kn Kt Gn Gt Mu Bn Bt Bm Eps",Kn,Kt,Gn,Gt,0.0,Bno,0.0,0.0,-0.1*Eps);
    B.Set(-3,"Kn Kt Gn Gt Mu Bn Bt Bm Eps",Kn,Kt,Gn,Gt,0.0,Bno,0.0,0.0,-0.1*Eps);
    B.Set(-4,"Kn Kt Gn Gt Mu Bn Bt Bm Eps",Kn,Kt,Gn,Gt,0.0,Bno,0.0,0.0,-0.1*Eps);
    B.Set(-5,"Kn Kt Gn Gt Mu Bn Bt Bm Eps",Kn,Kt,Gn,Gt,0.0,Bno,0.0,0.0,-0.1*Eps);
    B.Set(-6,"Kn Kt Gn Gt Mu Bn Bt Bm Eps",Kn,Kt,Gn,Gt,0.0,Bno,0.0,0.0,-0.1*Eps);
    B.Set(-7,"Kn Kt Gn Gt Mu Bn Bt Bm Eps",Kn,Kt,Gn,Gt,0.0,Bno,0.0,0.0,-0.1*Eps);
    dom.SetProps(B);

    for (size_t i=0;i<dom.BInteractons.Size();i++)
    {
        DEM::BInteracton * BI = dom.BInteractons[i];
        if (BI->P1==ParCtr||BI->P2==ParCtr) BI->valid = false;
    }

    String fkey_a(filekey+"_a");
    String fkey_b(filekey+"_b");

    // Here are some important variables to calculate the pressure in time
    dat.V0  = (dat.px1->x(0)-dat.px0->x(0)-dat.px1->Props.R+dat.px0->Props.R)*
              (dat.py1->x(1)-dat.py0->x(1)-dat.py1->Props.R+dat.py0->Props.R)*
              (dat.pz1->x(2)-dat.pz0->x(2)-dat.pz1->Props.R+dat.pz0->Props.R);
    dat.V0f = ParCtr->Props.V;
    dat.Vf  = ParCtr->Props.V;
    dat.ev  = 0.0;
    
    dom.Solve     (/*tf*/To   , /*dt*/dt, /*dtOut*/dtOut*10.0, &Setup1, &Report, fkey_a.CStr(),2,Nproc);

   
    dat.ev0 = dat.ev;
    dat.Vref= (dat.px1->x(0)-dat.px0->x(0)-dat.px1->Props.R+dat.px0->Props.R)*
              (dat.py1->x(1)-dat.py0->x(1)-dat.py1->Props.R+dat.py0->Props.R)*
              (dat.pz1->x(2)-dat.pz0->x(2)-dat.pz1->Props.R+dat.pz0->Props.R);
    // Here are some important variables to calculate the pressure in time
    dat.V0f = ParCtr->Props.V;
    dat.Vf  = ParCtr->Props.V;
    dat.V0  = dat.Vref;


    
    dom.Solve     (/*tf*/Tf+To, /*dt*/dt, /*dtOut*/dtOut, &Setup2, &Report, fkey_b.CStr(),2,Nproc);

}   
MECHSYS_CATCH


